Computer simulation of fatigue under diametrical compression 
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We study the fatigue fracture of disordered materials by means of computer simulations of a 
discrete element model. We extend a two-dimensional fracture model to capture the microscopic 
mechanisms relevant for fatigue, and we simulate the diametric compression of a disc shape speci- 
men under a constant external force. The model allows to follow the development of the fracture 
process on the macro- and micro- level varying the relative influence of the mechanisms of damage 
accumulation over the load history and healing of microcracks. As a specific example we consider 
recent experimental results on the fatigue fracture of asphalt. Our numerical simulations show that 
for intermediate applied loads the lifetime of the specimen presents a power law behavior. Under the 
effect of healing, more prominent for small loads compared to the tensile strength of the material, 
the lifetime of the sample increases and a fatigue limit emerges below which no macroscopic failure 
occurs. The numerical results are in a good qualitative agreement with the experimental findings. 

PACS numbers: 46.50.+a,62.20.Mk,05.10.-a 



I. INTRODUCTION 

Understanding the fatigue damage process of disor- 
dered materials is obviously important in many areas. 
The time scales related to the failure process, the rel- 
evant damage mechanisms, and scaling laws associated 
with this phenomenon have been widely studied both 
experimentally and theoretically p], 0, 01, 0]. Universal 
laws have been found in the fracture process related to 
both system sizes and geometry [5|, [6| , and in particular 
damage accumulation has been shown to play an impor- 
tant role in the formation and development of cracks with 
fractal structure [7|. 

This phenomenon is also important in practical ap- 
plications, in particular when dealing with asphalt mix- 
tures subject to traffic loading. Fatigue cracking is one 
of the main causes for asphalt layer failure in pavement 
structures, among moisture damage and thermal crack- 
ing. Most of the knowledge about fatigue failure of as- 
phalt mixtures relies, however, on experimental observa- 
tions [1, d, E3, EH, E2, EH • Material modeling to improve 
structural design is becoming an important tool to over- 
come this mode of material distress. 

Recently, fatigue life tests have been carried out to 
study the performance of asphalt mixtures [lj], by mea- 
suring the accumulation of deformation with time for 
cylindrical discs under diametrical compression applied 
periodically with a constant amplitude. Experiments re- 
vealed that the fatigue process has three different regimes 
depending on the external load amplitude a: when a 
falls close to the tensile strength of the material a rapid 
failure occurs, while for low load values a so-called fa- 
tigue limit emerges below which the material does not 
suffer macroscopic breaking. In the intermediate load 
regime the lifetime of the specimen tf has a power law 
dependence on the load, also called the Basquin law of 
fatigue [l5|]. Different models have been developed to 



obtain a theoretical understanding of the fatigue perfor- 
mance of asphalt specimens jig, Ell, Us). Besides dam- 
age accumulation, experimental research has shown the 
importance of the hea ling process to the lifetime of as- 
phalt mixtures @, E3, EI[ El, EH • The healing mecha- 
nism is related to the recovery of microcracks due to the 
viscoelastic nature of the binder material in polymeric 
mixtures, resulting in an extended lifetime. Theoreti- 
cal approaches turned out to have difficulties in captur- 
ing all the mechanisms involved in the fatigue process, 
such as damage growth, relaxation due to viscoelasticity 
and healing Q3, [IS EE E3, EI|. Although, these mod- 
els agree well with the results of the experiments, taking 
into account the stochastic nature of the fracture process 
and the cumulative effect of the load history, they fail 
to describe explicitly the processes of damage, fracture 
and failure of the solid. A more complete understand- 
ing of the failure process due to damage accumulation in 
the Brazilian test configuration may benefit from detailed 
computer simulations confronted with the experiments, 
as already was successfu lly per formed in many fracture 
processes @, [H, [H, EE MM HE M 13 • 

In this paper we enhance a 2D discrete element model 
of the fracture of disordered materials in order to cap- 
ture the relevant microscopic mechanisms of fatigue. In 
the model convex polygons symbolizing grains of the ma- 
terial are coupled by breakable beams, which can suffer 
stretching and bending deformation 0, 0, [H, 0, EE ED, 
l32l l33l l34i l35l |36| . The beams of the model can repre- 
sent the polymeric binder of asphalt between aggregates. 
Similarly to former studies, breaking of a beam can be 
caused by stretching and bending captured by the von 
Mises-type breaking criterion [22LI24I [25l [35[ . Addition- 
ally, we introduce an ageing mechanism, i.e. intact beams 
undergo a damage accumulation process which can again 
give rise to breaking. The accumulation process intro- 
duces memory over the loading history of the system. 
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FIG. 1: (Color online) Geometrical layout of the model. A 
disk shaped sample is subjected to diametrical compression. 
The disk is composed of convex polygons connected by elastic 
beams (magnified image). 



Damage recovery leading to healing is implemented in 
the model by limiting the range of memory which con- 
tributes to the amount of damage of fibers. We study in 
details the time evolution of the system and demonstrate 
that the model provides a good quantitative description 
of the experimental findings, where damage accumulation 
and healing proved to be essential. 



II. MODEL 

Our model is an extension of a realistic discrete ele- 
ment model of disordered materials which has been suc- 
cessfully applied to study various aspects of fracture and 
fragmentation phenomena [22|, 0, EjJ [35[ . The material 
is assumed to be composed of a large number of meso- 
scopic elements that interact elastically with each other. 
In our model the mesoscopic elements are convex poly- 
gon s randomly generated using a Voronoi construction 

The mesoscopic elements are elastic and unbreakable 
bodies, all with the same physical properties. Defor- 
mations are represented in the model by the possible 
overlap of the elements when pressed against each other. 
Between overlapping polygons we introduce a repulsive 
force which is proportional to the element's bulk Young 
modulus Y and the overlapping ar ea, and has a direction 
perpendicular to the contact line [25|, [37[ . 

Cohesion between the elements is introduced in the 
model by the inclusion of beams between neighboring 
polygons, see Fig. dJ The use of beams, as opposed to 
simple springs, takes into account tension and bending 
deformations, which are especially important in the com- 



TABLE I: Parameter values used 


in the simulations. 


Parameter 




Value 


Number of elements 




5070 


Density 


P 


r / 3 

5 g/cm 


Bulk Young modulus 


Y 


1 x 10 iU dyn/cm z 


Beam Young modulus 


Hi 


5 x 10 u dyn/cm z 


Time step 


St 


1 v in -6 c 


Diameter of the disk 


D 


20 cm 


Typical size of a single element 




0.5 cm 


Width of the load platen 


b 


2.5 cm 


Memory factor 


fo 


10 - 500 


Range of memory 


T 


500 - oo 



plex multi-axis stress field regime present in the geome- 
tries used experimentally. This model also resolves the 
local torques and shear forces which naturally arise when 
dealing with materials with complex micro- structures, 
where grains or aggregates are embedded in a matrix 
material. Beams have been particularly used in the lit- 
erature to prevent an unphysical failure of the system 
under shear [35[ , which happens when using central force 
elements. A detailed description of the beam model is 
given elsewhere [25|, [3l|, [37| . Briefly, the centers of mass of 
neighboring elements are connected by beams which can 
suffer stretching and bending deformation, and hence, 
exert an elastic force and torque on the polygons when 
stretched or bent 0, H S 0, M OH 1MB S S • 

In the model beams can break according to specified 
rules in order to explicitly model damage, fracture and 
failure of the solid [25|, [37[ . The imposed breaking rule 
must take into account both the immediate breaking of 
the beams by stretching and bending, as well as the ac- 
cumulation of damage and the healing mechanism. In 
order to achieve these conditions for each beam, at time 
t, we evaluate the quantities p(t) and q(t) defined by 

p(t)=f-g-V+ max( ' gi| ' |gjl) , (la) 

q(t) = p(t) + f f e Z ^ n p(t')dt f , (lb) 
Jo 

where e = Al/lo is the longitudinal deformation of the 
beam, and 6i and 0j are the rotation angles at the ends 
of beams between sites i and j, respectively. Equation 
([Taj) has the form of the von Mises plasticity criterion de- 
scribing the mechanical strength of beams with respect to 
stretching and bending deformation. The first part of Eq. 
([Taj) refers to the breaking of the beam through stretching 
and the second through bending, with e t h and 6 t h being 
the threshold values for elongation and bending, respec- 
tively. Equation (|lb|) refers to the long term memory, 
and accounts for the accumulation of damage and the 
healing mechanism. The rate of damage accumulation is 
assumed to have the form Ac(t) = fop(t) : i.e. the dam- 
age accumulated until time t is obtained as an integral 
over the loading history of elements c(t) = fo f*p(t')dt' . 
In our model healing is captured such that microcracks 



nucleated in beams can recover, a process which limits 
the total amount of damage. We describe this effect by 
introducing a finite range of memory over which the load 
experienced by beams has a contribution to the dam- 
age. For polymeric materials the dynamics of long chain 
molecules typically lead to an exponential form of the 
healing term. In Eq. (|lb|) the parameter fo is a memory 
factor and r is the time range of the memory over which 
the loading history of the specimen contributes to the ac- 
cumulation of damage 0] , and therefore a parameter that 
controls the healing mechanism. The breaking condition 
Eqs. (fTa|) and (fTbj) are evaluated at each iteration time 
step and those beams for which q(t) > 1 are considered 
to be broken, i.e., are removed from the simulation. 



For simplicity, all the beams have the same thresh- 
old values Sth an d Oth, and disorder is introduced in the 
model solely through the mesh generation [22|, [23|, 0, 
l25l l37l [38[ . The global material properties can be tuned 
by adjusting geometrical and microscopic physical pa- 
rameters of the model. The randomness of the Voronoi 
tessellation and the average size of the mesoscopic ele- 
ments, which define the average length and width of the 
beams are geometrical parameters that, together with 
the Young modulus of polygons Y and that of the beams 
E determine the macroscopic response of the model ma- 
terial 0, [H, [H, El, Ezl S|. The time evolution of the 
system is calculated by numerically solving the equations 
of motion for the translation and rotation of all elements 
using a sixth-order predictor-corrector algorithm. The 
breaking rule is evaluated at each iteration step. The 
breaking of beams is irreversible, which means that it is 
excluded from the force calculations for all following time 
steps. Table U summarizes the parameter values used in 
the simulations. 



Since usually the fracturing of materials is highly de- 
pendent on the structural environment in which it is stud- 
ied, one needs to reproduce the experimental geometry 
as accurately as possible. Starting from the Voronoi con- 
struction of a square, one cuts a disk-shaped specimen 
by removing those polygons outside the disk region and 
reshaping those at the border. The loading conditions 
for the numerical specimen are shown in Fig. [TJ This 
specimen is loaded in diametral compression by non- 
deformable platens, which are modeled by introducing 
two new polygons with the same physical properties of 
the loaded material. The load platens are cut with the 
same curvature of the disk. The ratio between the width 
of the load platen b and the diameter of the disk D is 
b/D = 0.125, as in the experiments of Ref. P3]. In the 
simulations presented in this paper each disk specimen is 
composed of 5070 mesoscopic elements. 
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FIG. 2: (Color online) Fracture of a disk-shaped solid sub- 
jected to a constant diametric load. Snapshots of the evolving 
damage. Simulations were carried out with the parameter 
values: a/a c = 0.8, f = 100, and r = 1300. For further 
parameters see Table |U 



III. RESULTS AND DISCUSSION 
A. Fracture development 

In the numerical experiments the sample is loaded with 
a constant external stress a. In order to avoid undesir- 
able large elastic waves, the load applied to the opposite 
platens is slowly increased from zero. When the load 
reaches the value a, the numerical sample is allowed to 
equilibrate during thousands of time steps including a 
small damping between the grains and friction according 
to Coulomb's friction law [25|, [37[. Only after equilibra- 
tion the breaking rules are applied to the beams. With 
this slow loading and long equilibration time, the vibra- 
tions of the sample are drastically reduced compared to 
the case when a constant load is applied instantaneously 
to the platens. The parameter values used in the simu- 
lations are summarized in Table H 

First we carried out computer simulations in order to 
determine the quasi-static strength a c of the model solid 
at the given parameters. Then the damage and fracture 
of the disc was studied varying the external load a below 
a c . In order to understand the failure process we studied 
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FIG. 3: ( Color online) Fraction of broken beams due to imme- 
diate breaking (squares) and memory accumulation (circles) 
as a function of a, for fo = 100, r = 500, together with the to- 
tal fraction of broken beams (diamonds). The insets show the 
final fracture pattern for a = 0.2 (left) and a = 0.8 (right). 




FIG. 4: (Color online) Final fracture pattern with fo = l (a) 
and fo = 100 (b) obtained by computer simulations at the 
same load a /a c = 0.8 and r — > oo. At the higher value of the 
memory factor fo more dispersed cracking occurs (6), while 
for low fo correlated crack growth can be obtained due to the 
dominance of immediate breaking. 



in details the evolution of the crack pattern during the 
simulations. Figure [2] presents snapshots of the evolving 
system for a disk of diameter 20 cm (approximately 80 
polygons) and a/a c = 0.8. In the initial stage of fracture, 
shown in Fig. [2](a), a number of cracks appears in the 
regions situated near the border of the platens contacts. 
Because of the finite width of the platens, the stress field 
induced in the sample is not purely uniaxial tensile stress 
along the load direction as one might expect [40[ , so the 
cracks do not initiate from the middle of the specimen. 
In fact, the cracks initiate from the border and coalesce 
to form wedges, as shown in Figure [2jb). The wedges 
extend up to approximately 0.25D, exhibiting a series of 
parallel vertical cracks, characteristic of compressive tests 
[251 ]. As time evolves the wedges penetrate the specimen 
(see Fig. Hfc)) originating a fracture in the center of the 
specimen, in a narrow region between the load platens. 
This mechanism leads to the sudden and catastrophic 
failure of the specimen (as shown in Figure [2] (d)). This 
failure behavior agrees qualitatively well with previous 
experimental observations [14|, [26[ . 

It is important to emphasize that without the time 
dependent damage accumulation, the system would not 
have macroscopic failure under a load a < a c . After 
some cracking events the disc would attain equilibrium 
resulting in an infinite lifetime. However, according to 
the failure criterion Eq. ([lb)) , intact elements accumulate 
damage which results in beam breaking after a finite time 
even under a constant local load. These beam breakings 
due to damage accumulation give rise to load redistri- 
bution and stress enhancements on intact elements trig- 
gering additional breakings. In this way, if healing is 
neglected, or analogously r — > oo is set in Eq. (|lb|h the 
system would fail after a finite time at any small load 
value. One can conclude from Eq. ([Tbj) that the main 



role of the memory factor /o is to control the relevant 
time scale of the system, i. e. the lifetime of the specimen 
tf has a dependence tf ~ 1 /Jo- 
in order to analyze the importance of the accumulation 
of damage in the fracture process, we have investigated 
the final fractions of the beams that have been broken 
during the course of the simulations due to the immediate 
breaking rule (first term in Eq. (|lbp) and those that have 
been broken due to damage accumulation (second term 
in Eq. (|lbp ) for different applied loads. Figure [3] shows 
the fractions of the broken beams due to the different 
breaking modes as a function of the applied load. For low 
load values damage dominates the failure of beams, while 
immediate breaking plays only a minor role. As the ex- 
ternal load is increased, the fraction of broken beams due 
to damage accumulation decreases monotonically, while 
that due to immediate breaking increases and saturates. 
For (7 > 0.65<j c the immediate breaking mode becomes 
dominant. It is interesting to note that the total fraction 
of broken elements is a decreasing function of the load, 
implying that correlated crack growth dominates when 

One of the most important features of our modeling ap- 
proach is that it accounts for the inhomogeneous stress 
field in the specimen and provides also information on 
the microstructure and spatial distribution of damage. 
The ageing of beams is less sensitive to the details of the 
stress field in the specimen. That is why when damage 
accumulation dominates the breaking of beams microc- 
racks are more dispersed in the stripe of the specimen 
between the loading platens. The insets of Fig. [3] show 
the final fracture pattern corresponding to a — 0.2a c and 
a = 0.8<j c . For a = 0.2 (left), where the effect of fatigue 
is more important, we observe more disordered cracks, 
while for a = 0.8, where the immediate breaking process 
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FIG. 5: Deformation e of the disk as a function of time t for 
fo = 100 and r = oo varying the external load a. Macroscopic 
failure is preceeded by an acceleration of the deformation. 
The vertical dashed line indicates the lifetime tf of the system. 
Note that tf decreases with increasing a. 
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FIG. 6: Fit of the experimental results of Ref. [lj]. Using 
the parameter values fo = 100 and r = 1300 a good quality 
agreement is obtained with the experiments. 



is more important, longer cracks can be identified due to 
crack growth. 

In our model the importance of damage accumulation 
for the beam breaking is controlled by the memory factor 
fo with respect to the immediate failure of elements. Fig- 
ure 2] demonstrates crack patterns obtained at the same 
load without healing r — > oo varying the value of Jo- 
lt can be observed that for large values of /o, even at 
a relatively high load a/a c = 0.85, a large amount of 
distributed cracking occurs in the specimen. Correlated 
crack growth can only be observed for low fo (see Fig.Uk) 
where immediate breaking has dominance. These crack 
patterns are very similar to those originated due to mem- 
ory effects in thermal fuse networks [4lj. The effect of 
increasing the memory factor fo resembles the effect of 
decreasing the applied load, in the sense that increasing 
fo has the same effect of enlarging the lifetime of the 
sample, as can be inferred from Eq. (|lb|) . 



B. Macroscopic time evolution 

On the macroscopic level the time evolution of the sys- 
tem is characterized by the overall deformation s of the 
disk and the lifetime tf of the specimen. The total de- 
formation of the specimen, £, is obtained directly from 
the simulation by monitoring the position of the loading 
platens. Figure [5] shows the behavior of e as a func- 
tion of time for different values of a with respect to the 
tensile strength <j c of the simulated specimen. Due to 
accumulation of damage, e increases monotonically until 
catastrophic failure of the material occurs. Some fluctu- 
ations can be seen in s before it starts increasing rapidly, 
due to the natural randomness of the beam breaking. 
The lifetime tf of the specimen is defined as the time 
at which the deformation e diverges. It can be observed 



that, by increasing a, the lifetime tf rapidly decreases 
and actually goes to zero as the external load approaches 
the tensile strength of the material <j c . 

For the purpose of fitting the experimental results, the 
most important parameters of the model that can be 
tuned are the two threshold values e t h and t h for the 
stretching and bending deformation of the beams, as well 
as the memory factor fo and the range of memory r. As 
we have pointed out, fo only sets the time scale of the 
system, while the other three parameters have a strong 
effect on the qualitative form of e(t) and tf(a). 

In the experiments of Ref. [14] the asphalt specimen 
was subjected to periodic loading with a constant ampli- 
tude monitoring the maximum deformation as a function 
of the number of loading cycles and the number of cycles 
to failure. In our model a constant load or a periodic 
load with a constant amplitude would have the same ef- 
fect, hence, we directly compare the simulation results to 
the experimental findings. 

Figure [6] shows a fit of the experimental results of Ref. 
[HI to the simulations, obtained by relating time to num- 
ber of cycles and an appropriate choice of geometrical 
and physical parameters for the mesoscopic elements and 
beams. We see that the simulation and the experimental 
curves show excellent agreement. 

In order to examine the consequences of the predom- 
inance of stretching or bending breaking modes we per- 
formed a number of simulations with different threshold 
values. The dependence of the lifetime on the external 
applied load is shown in Figure \J\ for r = oo and dif- 
ferent values of e t h and t h- The simulations show that 
for external loads a — » <r c , for all curves the lifetime de- 
creases rapidly and approaches the immediate failure of 
the specimen. For intermediate values of applied load, 
tf exhibits a power law behavior, in accordance with the 
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FIG. 7: Lifetime as a function of the load a for different 
breaking thresholds. In (a) the stretching threshold is kept 
fixed Eth — 0.03 and we vary the bending angle threshold 
9th, while in (b) th — 15° and the lifetime is calculated for 
different values of e t h • 



Basquin law of fatigue [if 



(2) 



where the value of the exponent 7 can be controlled by 
the breaking parameters Sth^th- We can see in Fig. [T^i 
that for a fixed value of eth — 0.03 the Basquin exponent 
varies from 7 = —1.1 ±0.1 for 0th = 3° (bending breaking 
mode dominance) to 7 = — 1.6 ± 0.1 as we increase 0th 
to 20°, hence increasing the influence of the stretching 
breaking mode. Obviously, the value of a c also increases 
with increasing 0th, since one of the breaking modes gets 
completely suppressed. Similarly in Fig. [7)3, with a fixed 
value of t h = 15° the Basquin exponent varies from 7 = 
-1.4 ± 0.1 for e t h = 0.04 to 7 = -1.8 ± 0.1 for e th = 
0.01, in which case there is a clear predominance of the 



FIG. 8: Lifetime as a function of the load a/a c . The square 
filled symbols correspond to experimental results, the open 
circles correspond to r = 00, and the *, + and x symbols to 
r = 21150, 16450 and 11750, respectively. 



stretching breaking mode. 

In order to obtain the measured value 7 = 2.0 ± 0.1 
the breaking parameters had to be set such that the 
stretching mode of Eq. ([Taj) dominates the breaking, i.e. 
eth — 0-01 an d t h = 20° were used which implies that 
the bending mode has only a minor role in breaking. The 
upturn of the curves of tf for small loads in Fig. [8] indi- 
cates the emergence of the fatigue limit 07 in the system 
below which a < 07 the specimen only suffers partial fail- 
ure and has an infinite lifetime. Figure [8] also illustrates 
that the value of o\ is determined by the range of memory 
r, since healing has only a dominating effect on the time 
evolution of the system when r is comparable to the life- 
time of the specimen measured without healing (r — > 00). 
Hence, it follows from the Basquin law, Eqj2j that the 
fatigue limit 07 scales with r as 07 ~ r~~ . The good 
quantitative agreement between simulations and experi- 
ment obtained for r = 21150 indicates that the stretching 
breaking mechanism is more important for this particular 
material. This is reasonable since the polymer binding 
in asphalt can not sustain torsion. 



IV. CONCLUSIONS 

We carried out a computational study of the fatigue 
failure of disordered materials occurring under a constant 
external load. A two-dimensional discrete element model 
was extended to capture microscopic failure mechanisms 
relevant for the process of fatigue. As a specific example, 
we considered experiments on asphalt where damage re- 
covery in the form of healing is known to play a crucial 
role for the long term performance of the material. 

The breaking of beams is caused by two mechanisms: 
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immediate breaking occurs when the failure thresholds 
of stretching and bending are exceeded. Intact beams 
undergo a damage accumulation process which is limited 
by the finite range of memory (healing). The relevant 
parameters of the model to fit the experimental results 
are the breaking thresholds of stretching e t h and bending 
t h and the range of memory r. We carried out a large 
amount of computer simulations to study the time evo- 
lution of the system at the macro- and micro-level. We 
demonstrated that at the micro- level failure of beams due 
to damage is responsible for the diffuse crack pattern, 
while correlated crack growth is obtained when immedi- 
ate breaking has dominance. The model provides a good 
quantitative agreement with previous experimental find- 
ings on the fatigue failure of asphalt, i.e. varying only 
three parameters eth, Oth, and r of the model we could re- 
produce the deformation-time diagram and the Basquin- 
law of asphalt obtained in experiments. Very interest- 
ingly we found that the value of the Basquin exponent of 



the model is controlled by the relative importance of the 
stretching and bending modes in beam breaking. The 
fatigue limit o\ is obtained as a threshold value of the 
external load below which the specimen does not suffer 
macroscopic failure and has an infinite lifetime. In the 
framework of our model g\ is solely determined by the 
range memory r over which local loads contribute to the 
total accumulated damage in the system. 
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